function [toframes, todescr]=flow_singular_pattern_detect_low_memory(uv,varargin)
%   function [toframes, todescr,gss,ess]=flow_singular_pattern_detect(uv,varargin)
%       SIFT-like flow field singular-pattern detector. 
%
%       See SIFT for the details of options.
%
%       OUTPUT:
%                   toframes - positions, scales, and singular energy of
%                                       the detected singular patterns.
%                               toframes(1:2,:) - spatial coordinates
%                               toframes(3,:)   -   scales (or radius) 
%                               toframes(4,:)   -   singular energy
%                   todescr    - descriptor vectors
%                   gss           - gaussian cascade of the flow field
%                                      (octaves)
%                   ess           - singular energy cascade
%       
%       by Wei Liu @ Florida Tech 2010
%                           Adapted from Andrea Vedaldi's implementation of
%                           SIFT detector
%
%
%
%
%
% This file is part of the source-code demonstrating the method published 
% in the paper: 
%
% Wei Liu and Eraldo Ribeiro. Scale and Rotation Invariant Detection 
%      of Singular Patterns in Vector Flow Fields. IAPR International 
%      Workshop on Structural Syntactic Pattern Recognition (S-SSPR), 
%      Cesne, Turkey, 2010.
%
% BIBTEX ENTRY:
%
% @inproceedings{WeiRibeiroSSPR2010,
%   author    = {Wei Liu and Eraldo Ribeiro},
%   title     = {Scale and Rotation Invariant Detection of Singular 
%                Patterns in Vector Flow Fields},
%   booktitle = {IAPR International Workshop on Structural Syntactic 
%                Pattern Recognition (S-SSPR)},
%   year      = {2010},
%   pages     = {XXX-XXX}
% }
%
% The source-code can be obtained from 
%        http://www.cs.fit.edu/~eribeiro/flowdetector/
% 
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you 
% maintain the name of the original authors.
%
% (c) 2010 by Wei Liu and Eraldo Ribeiro 
% eribeiro@cs.fit.edu
% Florida Instititute of Technology, Melbourne, Florida, U.S.A.
%







%%
if(nargin<1)
    error('At least one argument is required.');
end
%%

verb=1;
[M,N,C] = size(uv) ;

% Lowe's equivalents choices
S      =  3 ;
%S      =  5;
omin   =-1 ;
O      = floor(log2(min(M,N)))-omin-2 ; % up to 8x8 images
sigma0 = 1.6*2^(1/S) ;                  % smooth lev. -1 at 1.6
sigman = 0.5 ;
%thresh = 0.04 / S / 2 ;
thresh=0;
magnif = 3.0 ;
k=2^(1/S);
norder=1;
lorder=norder;
if verb>0, fprintf('SIFT: computing scale space...\n') ; tic ; end

% Parese input
discard_boundary_points = 0 ;
verb = 0 ;

for vv=1:2:length(varargin)
    switch lower(varargin{vv})
        case 'norder'                                           % maximum positive order
            norder=varargin{vv+1};
        case 'lorder'
            lorder=varargin{vv+1};
            
        case 'numoctaves'
            O = varargin{vv+1} ;
            
        case 'firstoctave'
            omin = varargin{vv+1} ;
            
        case 'numlevels'
            S = varargin{vv+1} ;
            
        case 'sigma0'
            sigma0 = varargin{vv+1} ;
            
        case 'sigman'
            sigman = varargin{vv+1} ;
            
        case 'threshold'
            thresh = varargin{vv+1} ;
            
        case 'boundarypoint'
            discard_boundary_points = varargin{vv+1} ;
            
        case 'magnif'
            magnif = varargin{vv+1} ;
            
        case 'verbosity'
            verb = varargin{vv+1} ;
                       
        otherwise
            error(['Unknown parameter ''' varargin{vv} '''.']) ;
    end
end

%%
if(norder<0||lorder<0)
    error('norder and lorder must both be positive');
end
%%
% this function takes very long time to return. Parallel implementation
% (e.g. GPU) would be good
if verb>0
    fprintf('computing gaussian cascade...lots of convolution\n');
    tic();
end

nfilters=(norder+lorder+1)*2;
% --------------------------------------------------------------------
%                                                  Check the arguments
% --------------------------------------------------------------------
I=uv;

smin=-1;smax=S+1;

toframes=[];todescr=[];
% --------------------------------------------------------------------
%                                                           Do the job
% --------------------------------------------------------------------

% Scale multiplicative step
k = 2^(1/S) ;

% Lowe's convention: the scale (o,s)=(0,-1) has standard deviation
% 1.6 (was it variance?)
if(nargin < 7)
    sigma0 = 1.6 * k ;
end

if omin < 0
    for o=1:-omin
        I = doubleSize(I) ;
    end
elseif omin > 0
    for o=1:omin
        I = halveSize(I) ;
    end
end

[M,N,nc] = size(I) ;

% Index offset
so = -smin+1 ;

% --------------------------------------------------------------------
%                                                         First octave
% --------------------------------------------------------------------
%
% The first level of the first octave has scale index (o,s) =
% (omin,smin) and scale coordinate
%
%    sigma(omin,smin) = sigma0 2^omin k^smin
%
% The input image I is at nominal scale sigman. Thus in order to get
% the first level of the pyramid we need to apply a smoothing of
%
%   sqrt( (sigma0 2^omin k^smin)^2 - sigman^2 ).
%
% As we have pre-scaled the image omin octaves (up or down,
% depending on the sign of omin), we need to correct this value
% by dividing by 2^omin, getting
%e
%   sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 )
%

if(sigma0 * 2^omin * k^smin < sigman)
    warning('The nominal smoothing exceeds the lowest level of the scale space.') ;
end

%SS.octave{1} = zeros(M,N,2,smax-smin+1) ;% first octave for a 2-channel image
ESS=zeros(M,N,smax-smin+1);
STACK_SS=zeros(M,N,nfilters,smax-smin+1);
fb=cell(smax-smin+1,1);
SS = zeros(M,N,2,1);
 SS(:,:,1,1)  = imsmooth(I(:,:,1), ...
     sqrt((sigma0*k^smin)^2  - (sigman/2^omin)^2)) ;
 SS(:,:,2,1)  = imsmooth(I(:,:,2), ...
     sqrt((sigma0*k^smin)^2  - (sigman/2^omin)^2)) ;%smooth to each channel

radius=round(sigma0*k^smin*4);
if verb>0
    fprintf('creating laurent basis flows ...\n');
    tic();
end

fb{1}=create_laurent_filters(radius,lorder,norder,sigma0*k^smin);
% In out hybrid approach, the scaling within an octave is applied to the
% basis function, not on the input image.
for s=smin+1:smax
    % Here we go from (omin,s-1) to (omin,s). The
    % standard deviation is
    %
    %  (sigma0 2^omin 2^(s/S) )^2
    % Aftred dividing by (2^omin)^2 (to take into account the fact
    % that the image has been pre-scaled omin octaves), the
    % standard deviation of the smoothing kernel is
    %
    %   dsigma = sigma0 k^s
    %
    dsigma = k^s * sigma0 ;
    radius=round(dsigma*4);
    %  fb{s+so}=create_monomial_filters(radius,norder,dsigma);
    fb{s+so}=create_laurent_filters(radius,lorder,norder,dsigma);
end

if verb>0
    fprintf('processing ocatve 1 \n');
    tic();
end

for s=smin:smax
    stack=apply_monomial_filters(SS,fb{s+so});
    dsigma = k^s * sigma0;
    fprintf('scale:%f\n',dsigma);
    stack=(stack/dsigma).*(2^(omin));%normalize the descriptor against the scale.
    STACK_SS(:,:,:,s+so)=stack;
    ESS(:,:,s+so)=total_singular_energy(stack,lorder);
end

idx = siftlocalmax(ESS, thresh) ;
[seng]=ESS(idx);%obtain singular energy

%     K=length(idx) ;
[i,j,si] = ind2sub( size( ESS ), idx ) ;
y=i ;
x=j ;
s=si-1+smin ;

%% we have to use this twisted code to avoid using PERMUTE(STACK_SS)
i=repmat(i,[nfilters,1]);
j=repmat(j,[nfilters,1]);
si=repmat(si,[nfilters,1]);
tidx=repmat([1:nfilters]',[1,size(i,2)]);
idx=sub2ind(size(STACK_SS),i,j,tidx,si);
descr=reshape(STACK_SS(idx),[nfilters, size(i,2)]);
oframes = [x(:)';y(:)';s(:)';seng(:)'] ;

if verb > 0
    fprintf('%d initial points (%.3f s)\n', ...
        size(oframes, 2), toc) ;
    tic ;
end

% Remove points too close to the boundary
if discard_boundary_points
    % radius = maginf * sigma * NBP / 2
    % sigma = sigma0 * 2^s/S
    
    %rad = magnif * gss.sigma0 * 2.^(oframes(3,:)/gss.S) * NBP / 2 ;
    rad=(magnif * sigma0 * 2.^(oframes(3,:)/S))/2;
    sel=find(...
        oframes(1,:)-rad >= 1                     & ...
        oframes(1,:)+rad <= size(GSS,2) & ...
        oframes(2,:)-rad >= 1                     & ...
        oframes(2,:)+rad <= size(GSS,1)     ) ;
    oframes=oframes(:,sel) ;
    descr=descr(:,sel);
    if verb > 0
        fprintf('SIFT: %d away from boundary\n', size(oframes,2)) ;
        tic ;
    end
end

%% calculate its position and radius
oframes(1:2,:)=oframes(1:2,:).*2^(omin);
dsigma=k.^oframes(3,:).*sigma0;
radius=dsigma*2^(1+omin-1)*1.5;
oframes(3,:)=radius;
toframes=cat(2,toframes,oframes);
todescr=cat(2,todescr,descr);

%%

% --------------------------------------------------------------------
%                                                        Other octaves
% --------------------------------------------------------------------

for o=2:O
    % We need to initialize the first level of octave (o,smin) from
    % (o-1,smin)
    %
    % The extra variance needed is:
    %  ( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(smin/S - 1) )^2
    %
    % As usual, we divide by 2^o to cancel out the effect of the
    % downsampling and we get
    %
    % sigma0 k^smin sqrt(1-1/k^2S)
    %
    if verb>0
        fprintf('processing octave:%d\n',o);
        tic();
    end
    
    dsigma=sigma0*(k^smin)*sqrt(1-1/(k^(2*S)));
    %% let's calculate the best matched flow field first
    
    %%
    TMP = halveSize(squeeze(SS(:,:,:,1))) ;
    [M,N,nc] = size(TMP) ;
    
    TMP(:,:,1)= imsmooth(squeeze(...
        TMP(:,:,1)),...
        dsigma);
    TMP(:,:,2)= imsmooth(squeeze(...
        TMP(:,:,2)),...
        dsigma);
    SS= zeros(M,N,2,1) ;
    SS(:,:,:,1) = TMP ;
    ESS=zeros(M,N,smax-smin+1);
    STACK_SS=zeros(M,N,nfilters,smax-smin+1);
    for s=smin:smax
        stack=apply_monomial_filters(SS,fb{s+so});
        dsigma = k^s * sigma0;
        stack=(stack/dsigma).*(2^((o+omin-1)));%normalize the descriptor against the scale.
        %stack=(stack/(dsigma));
        STACK_SS(:,:,:,s+so)=stack;
        ESS(:,:,s+so)=total_singular_energy(stack,lorder);
    end
    
    %% extract descriptors
    idx = siftlocalmax(ESS, thresh) ;
    [seng]=ESS(idx);%obtain singular energy
    
    %     K=length(idx) ;
    [i,j,si] = ind2sub( size( ESS ), idx ) ;
    y=i ;
    x=j ;
    s=si-1+smin ;
    
    %% we have to use this twisted code to avoid using PERMUTE(STACK_SS)
    i=repmat(i,[nfilters,1]);
    j=repmat(j,[nfilters,1]);
    si=repmat(si,[nfilters,1]);
    tidx=repmat([1:nfilters]',[1,size(i,2)]);
    idx=sub2ind(size(STACK_SS),i,j,tidx,si);
    descr=reshape(STACK_SS(idx),[nfilters, size(i,2)]);
    
    oframes = [x(:)';y(:)';s(:)';seng(:)'] ;
    
    if verb > 0
        fprintf('SIFT: %d initial points (%.3f s)\n', ...
            size(oframes, 2), toc) ;
        tic ;
    end
    
    % Remove points too close to the boundary
    if discard_boundary_points
        % radius = maginf * sigma * NBP / 2
        % sigma = sigma0 * 2^s/S
        
        %rad = magnif * gss.sigma0 * 2.^(oframes(3,:)/gss.S) * NBP / 2 ;
        rad=(magnif * sigma0 * 2.^(oframes(3,:)/S))/2;
        sel=find(...
            oframes(1,:)-rad >= 1                     & ...
            oframes(1,:)+rad <= size(GSS,2) & ...
            oframes(2,:)-rad >= 1                     & ...
            oframes(2,:)+rad <= size(GSS,1)     ) ;
        oframes=oframes(:,sel) ;
        descr=descr(:,sel);
        if verb > 0
            fprintf('SIFT: %d away from boundary\n', size(oframes,2)) ;
            tic ;
        end
    end
    
    %% calculate its position and radius
    oframes(1:2,:)=oframes(1:2,:).*2^(o+omin-1);
    dsigma=k.^oframes(3,:).*sigma0;
    radius=dsigma*2^(o+omin-1)*1.5;
    oframes(3,:)=radius;
    toframes=cat(2,toframes,oframes);
    todescr=cat(2,todescr,descr);
end

if(size(toframes,2)~=size(todescr,2))
    error('odd things happened');
end

% -------------------------------------------------------------------------
%                                                       Auxiliary functions
% -------------------------------------------------------------------------

function J=flowResize(I, scale)
   u=imresize(I(:,:,1),scale,'bilinear');
   v=imresize(I(:,:,2),scale,'bilinear');
   J=cat(3,u,v)*scale;

function J = doubleSize(I)
J=zeros(size(I)*2);
J(:,:,1)=doubleSizeHelper(I(:,:,1));
J(:,:,2)=doubleSizeHelper(I(:,:,2));
J=J*2;

function J=doubleSizeHelper(I)
[M,N]=size(I) ;
J = zeros(2*M,2*N) ;
J(1:2:end,1:2:end) = I ;
J(2:2:end-1,2:2:end-1) = ...
    0.25*I(1:end-1,1:end-1) + ...
    0.25*I(2:end,1:end-1) + ...
    0.25*I(1:end-1,2:end) + ...
    0.25*I(2:end,2:end) ;
J(2:2:end-1,1:2:end) = ...
    0.5*I(1:end-1,:) + ...
    0.5*I(2:end,:) ;
J(1:2:end,2:2:end-1) = ...
    0.5*I(:,1:end-1) + ...
    0.5*I(:,2:end) ;

function J = halveSize(I)
J=I(1:2:end,1:2:end,:) ;
J=J*0.5;

function [J]=total_singular_energy(stack,lorder)
stack=permute(stack,[3 1 2]);
J=sum(stack((lorder+1)*2+1:end,:,:).^2,1);
if(lorder>0)
    J=J+sum(stack(1:lorder*2,:,:).^2,1);
end
return;

